function [x, fx, exitflag, output] = bisezione(fun, a, b, options, varargin)
    % sepzeri   Ricerca sottointervalli con variazione di segno di f(x)
    %
    % Sintassi: Xrad
    % Impostazioni di default
    if nargin < 4
        xtol = 1e-8; ftol = 1e-6; maxiter = 100;
        verbose = false;
    else
        xtol = options.xtol; ftol = options.ftol;
        maxiter = options.maxiter; verbose = options.verbose;
    end
    
    % Inizializzazione variabili d'uscita
    x = []; fx = []; exitflag = -1;
    output = struct('iter', [], 'a', [], 'b', [], 'xm', [], 'fxm', []);
    
    % Verifica degli estremi dell'intervallo
    if (a >= b)
        error(sprintf('Metodo di bisezione non applicabile (%f>=%f)!', a, b));
    end
    
    % Re-imposta le tolleranze se troppo piccole
    xtol = max(xtol, 5*eps);
    ftol = max(ftol, 5*eps);
    
    % Funzione agli estremi
    fa = feval(fun, a, varargin{:});
    fb = feval(fun, b, varargin{:});
    
    % Verifia della funzione agli estremi
    if (fa>0) == (fb>0)
        error(sprintf('Metodo di bisezione non applicabile (fa * fb > 0) in [%f, %f]!', a, b));
    end
    
    % Costanti per il test di convergenza
    xref = abs(b-a);
    fref = max(abs([fa, fb]));
    
    % Inizializzazioni
    iter = 0; % Numero corrente di iterazioni
    
    % Output viene valutata solo se necessario
    if nargout == 4
        output.iter = 0;
        output.a = repmat(NaN, 1, maxiter);
        output.b = repmat(NaN, 1, maxiter);
        output.xm = repmat(NaN, 1, maxiter);
        output.fxm = repmat(NaN, 1, maxiter);
    end
    
    % Stampa informazioni se disponibili
    if verbose
        fprintf('%30s\n%4s %12s %12s\n%30s\n', ...
            repmat('-', 1, 30), 'iter', 'xm', 'f(xm)', repmat('-', 1, 30));
    end
    
    % Ciclo principale di convergenza
    while iter <= maxiter
        % Aggiorno contatore delle iterazioni
        iter = iter + 1;
        
        % Calcolo del punto medio
        dx = (b-a)/2;
        xm = a + dx;
        
        % Valuto la funzione nel punto medio
        fxm = feval(fun, xm, varargin{:});
        
        % Salvataggio informazioni se richieste
        if nargout == 4
            output.iter = iter
            output.a(iter) = a;
            output.b(iter) = b;
            output.xm(iter) = xm;
            output.fxm(iter) = fxm;
        end
        
        % Stampa informazioni se disponibili
        if verbose
            fprintf('%4d %12.4e %12.4e \n', iter, xm, fxm);
        end
        
        % Test di convergenza
        if (abs(fxm)/fref < ftol) || (abs(dx)/xref < xtol)
            % Trovata radice ed esci dal ciclo
            x = xm;
            fx = fxm;
            exitflag = 1;
            break;
        end
        
        % Aggiorno gli estremi dell'intervallo in base al segno della
        % funzione
        if sign(fa) == sign(fxm)
            % Zero tra [xm, b]
            a = xm;
            fa = fxm;
        else
            % Zero tra [a, xm]
            b = xm;
            fb = fxm;
        end
    end
    
    % Sistemazione dell'output finale
    if nargout == 4
        output.a = output.a(1:output.iter);
        output.b = output.b(1:output.iter);
        output.xm = output.xm(1:output.iter);
        output.fxm = output.fxm(1:output.iter);
    end
    
    if verbose
        fprintf('%30s\n', repmat('-', 1, 30));
    end
end

